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SUMMARY 


An exact, full- potential- equation (FPE) model for the steady, irrotational, ho- 
mentropic and homoenergetic flow of a compressible, homocompositional, inviscid 
fluid through two-dimensional planar cascades of airfoils was derived, together with 
its appropriate boundary conditions. 

A computer program, CAS2D, was developed that numerically solves an arti- 
ficially time- dependent form of the actual FPE. The governing equation was dis- 
cretized by using type- dependent, rotated finite differencing and the finite area tech- 
nique. The flow field was discretized by providing a boundary- fitted, nonuniform 
computational mesh. The mesh was generated by using a sequence of conformal 
mapping, nonorthogonal coordinate stretching, and local, isoparametric, bilinear 
mapping functions. The discretized form of the FPE was solved iteratively by using 
successive line overrelaxation. The possible isentropic shocks were correctly cap- 
tured by adding explicitly an artificial viscosity in a conservative form. In addition, 
a four- level, consecutive, mesh refinement feature makes CAS2D a reliable and fast 
algorithm for the analysis of transonic, two-dimensional cascade flows. 


INTRODUCTION 

This work is an extension of the author's doctoral research (ref. 1). 

The purpose of this report is to give instructions to potential users of the com- 
puter program CAS2D. These instructions refer to the possible applications and re- 
strictions of CAS2D. The basic assumptions of the theory used to develop CAS2D are 
also detailed. 

The simplest (although, still exact) mathematical model of the fluid flow is the full 
potential equation (FPE). Expressed in the two-dimensional Cartesian coordinate 
system (fig. 1), the FPE can be written as 

f 2 - XX + (a 2 - yy ’ 2 * x^, y* xy = 0 «> 

where a is the local speed of sound and cp is the velocity potential function (V = V<p). 
This equation represents a continuity equation for the steady, irrotational, homen- 
tropic flow of an inviscid, homocompositional, compressible fluid. The conditions for 


its validity are that no sources of entropy production exist and that the flow is adia- 
batic. This can be seen from Crocco* s formula 

V x (V x V) = VH - T VS (2) 

where V = ue x + v£y is the velocity vector; H, the total enthalpy; T, the absolute 
temperature; S, the entropy; and 


V x V= 0 


( 3 ) 


the irrotationality condition. 

In practice, this means that the flow is uniform (irrotational) at upstream infinity 
and that all the viscous effects are confined to a narrow, nonseparated boundary- layer 
region. The flow at downstream infinity should also be uniform in accordance with 
the steady- flow inviscid theory. 

In actual viscous flows the vorticity created at the solid boundary is convected and 
diffused downstream, where it forms a wake. The velocity deficit in the wake rapidly 
disappears. Therefore it is reasonable to assume that the downstream infinity is only 
a few chord lengths from the trailing edge. 

The discontinuities possible in the solution of FPE are isentropic shocks. They 
do not represent physical shock waves (ref. 2) because they do not satisfy Rankine- 
Hugoniot jump conditions. Nevertheless they are of approximately correct strength 
and at approximately the correct position if the Mach number just ahead of the dis- 
continuity is less than about 1. 3. 

The results obtained by this two-dimensional analysis are not directly applicable 
to three-dimensional, potential, rotating flows through a cascade of blades. The rea- 
son is that the Coriolis force does not exist in two-dimensional, potential, cascade 
flows. 


ANALYSIS 

Transonic, Full Potential Equation 

The continuity equation for steady flow with no sources or sinks in the flow field is 


V* (pV) = V • (Vp) +p(V • V) = 0 

where p is the local fluid density and V is the local fluid velocity vector 


( 4 ) 
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For continuous, homentropic flow of a homocompositional fluid the relation 


Vh _ Vp 

a 2_ P 


( 6 ) 


is valid, where h is the local static enthalpy and a is the local speed of sound 


a 2 = 


/3p\ 

\9p Jt 


( 7 ) 


isentropic 


After premultiplying equation (4) by a /p and taking into account equation (6), we get 


a 2 (V • V) + V • (Vh) = 0 


( 8 ) 


In general 


h = 



( 9 ) 


where the total enthalpy H must be constant throughout the flow field in order to pre- 
serve the irrotationality condition (eq. (2)). 

Then equation (8) becomes 


a 2 (V • V) - (V • V)Q = 0 


( 10 ) 


where 


Q = 


JvP 

2 


(ID 


The continuity equation (10) can be written in its quasi- conservative scalar form 
(ref. 3) as 


or in its fully conservative scalar form as 


_e_ 


a 




= 0 


(13) 


With introduction of the velocity potential, equation (10) becomes the full- potential 
equation 


(f 2 - xx + f 2 - <p*y )*, yy - 2 ? x* y* xy = 0 < 14 > 

Equation (14) is a quasi- linear, second-order partial differential equation of the 
mixed (elliptic-hyperbolic) type. Its canonical form (ref. 4) is 

(a 2 - qV H ss + a 2 (v 2 <p E - cp E gs ) = 0 (15) 

where the superscript H designates upstream differencing and the superscript E 
designates the central differencing that should be used in the discretization of equa- 
tion (15). Here, s is the streamline direction and 


q 2 = (u 2 + v 2 ) 


(16) 


For the finite difference evaluation of the derivatives, the flow field and the gov- 
erning equations were transformed into a rectangular (X, Y) computational space 
(fig. 4). The transformation matrix is 


[J] T = 


x , x y.x 


l x ,y y.Tj 


( 17 ) 


Then the modified contravariant velocity components are 



where D = det[J] and 


[B] =[JT 1 [j T ] 


-1 


(19) 


Consequently the fully conservative form of the continuity equation (eq. (13)) becomes 


((pU) j x + (pV) f yj ~ 0 


( 20 ) 


where 


P = 


^1+XZ-l M^d- q 2 )^ 


i/ (y- 1) 


and 


a 2 = 


i ( 1 -q 2 

WJ 2 


The quasi- conservative form (ref. 3) becomes 

a 2 (U x + V y ) - (UQ x + VQ Y ) = 0 (21) 

The principal part of the full- potential form (ref. 4) of the continuity equation 
(eq. (14) or (15)) transforms into 

E E E / jj Ei /H Ei 

S XX (p , XX + S YY < ? YY + S XY cp , XY + R XX XX _ (p , XX j + R YY YY " (p , Yy) 

+ r xy (^fxy ' ^fxy) = 0 (22) 

where the superscript H designates upstream (or backward) differencing and the 
superscript E designates central differencing. The coefficients of this equation are 
as follows: 
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u = max |l ; 0.1 

2 

q 


o _ „2 rtT} U 

S XX ” DB 11 ~ 


S YY - a DB 22 ^ 


S XY= 2 ( a2DE 12--^) 


>- 


(23) 


_ mU 

xx~ “d 


Ryy 


, r 2 

= _£iZ_ 




For a locally subsonic flow all the second derivatives will be evaluated by central 
differencing. Upstream differencing (as indicated in eqs. (15) and (22)) is a numer- 
ical attempt to simulate the exact shape of the upstream- facing domain of depend- 
ence of the locally hyperbolic full-potential equation. It can be shown (ref. 5) that 
this procedure introduces a numerical truncation error called artificial viscosity, 
which is proportional to 

j bSb 

This artificial viscosity term was added to equation (20) or (21) in a conservative 
form 


RVISC - f H x + g 1 ^ 


(24) 


where 
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(25) 


f i, j - ~^r ( 2R XX (? fxX + E XY ?i , E Xy) 

«i, j = - if (^fxY + 2 R YY*?Yy) < 26 > 

The conservative form of RVISC assures the unique position and strength (ref. 2) of the 
possible isentropic discontinuities in the solution of the governing equation. 

Boundary Conditions 

The boundary conditions for equation (14) (or (15)) that are valid in the present 
analysis of a two-dimensional cascade flow are the following: On the airfoil surface 
the condition is 


Vcp • VF = 0 


(27) 


where F = F(x, y) is the equation of the airfoil contour line. At axial infinities (x = 

±») the flow is assumed to be uniform. The upstream boundary condition is incor- 
porated in the special form of the potential function (ref. 4). 

cp(x, y) = (x cos o!_ 0O + y sin a _,*,) + G(x, y) (28) 

where a_ QO and are the angles between the free stream and the x axis at the 

upstream and downstream infinities, respectively. The term G(x,y) is a so-called 
M reduced potential” describing the disturbances with respect to the free stream at 
upstream infinity. All the velocities have been normalized with respect to the mag- 
nitude of the free- stream velocity vector at upstream infinity. Hence the boundary 


condition at upstream infinity (x = -°°) is 


(<P , X ) =COSQ! -oo 

(29) 

V ’ 1-00 


y)_„ = sln 

(30) 


That is, 
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(G) = Constant 

' — OO 


( 31 ) 


where the arbitrary constant is taken as zero. 

Conditions that should be applied along the identically shaped and periodically 
spatially positioned boundaries (corresponding to the lines cd and ab in fig. 1) 
represent the physical fact that the velocity vector is a periodic function in a cascade 
flow 


V(x, y) = V(x, y + h) 


According to the definition of the velocity potential 


cf> d = cp sl + f V<p • dT 
•'a 


where 



dy e. 


and because 



it follows that 


^d^a + f^ &<P 



or 


(<P d - <P a ) = h sin o!_ 0O 

From equation (18) it follows that 


(32) 


(33) 


(34) 


(35) 


(36) 


(37) 


<^d - <Ps) = < x d “ x a> cos “-«> + (yd - y a > sin a -oo + ( G d - G a > < 38 > 
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Because (x^ - x a ) = 0 and (y^ - y a ) = h, it follows from equations (27) and (28) that 
the boundary conditions along the periodic boundaries are 

G(x,y)= G(x,y +h) (39) 

The global mass conservation is expressed as (fig. 2) 

h.ooP.oo'Loo cos Q!_ 00 = cos (40) 

where p is density and q is the magnitude of the velocity vector. Because p_ oo = 1 
and q wOQ = 1, 

f(q^) = cos “-tco - cos ^-co = 0 ( 41 ) 


if h.^ = h^. Hence 


= P +oo COS “+00 + ^ cos «+» 


\*w 


(42) 


By using a simple Newton- Raphson iterative technique 


q n+ ^ = q* 

H -foo M +oo 


n 40 

f '(0 


(43) 


the value of q^ can be easily determined. Then the boundary condition at the down- 
stream infinity becomes 



COS a ^ - cos 


(44) 


For lifting flows the potential function becomes multivalued. This problem can be 
resolved by inserting an arbitrarily shaped cut of zero thickness in the flow field in 
such a way that it connects the airfoil with the infinity (fig. 1). This cut is conven- 
iently assumed to leave the trailing edge. The finite discontinuity in the velocity po- 
tential A <p at the trailing edge is equal to the circulation r of the velocity field. 

This constant discontinuity should be preserved (i. e. , there should be no jump in the 
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static pressure) at every point of the cut all the way to the downstream infinity. The 
exact value of r can be calculated from the input data and the value of q^ (eq. (33)). 

r = Acp = h(sin 0 !_ 00 - q^ sin (45) 

NUMERICAL METHOD 
Computational Mesh Generation 

The computational mesh (fig. 3) was generated by using a conformal mapping 
function of the form 


w = e 1 * 3 In + e-V In (46) 

\m + z/ \1 + mz / 

It conformally maps a unit circle with a slit in the middle (ref. 6) whose endpoints are 
situated at z = ±m onto the cascade of slits in the w = x + iy plane. 

The slits are spaced 2ir cos /3 apart, where /3 is the stagger angle of the cascade 
of slits (fig. 4). The unit circle with a slit in the z = | + irj plane can be M unwrapped 1 * 
by using elliptic polar coordinates (ref. 7): 


£ = m cosh u 0 cos v 

6 c 

(47) 

rj = m sinh u sin v 

6 G 

(48) 


The resulting (u , v ) plane can be transformed to the rectangular, computational (X, Y) 

6 G 

domain with the help of a nonorthogonal shearing transformation (fig. 4) of the general 
form 


Y = 


T(v e ) 


(49) 


X = v g + Y * DRD (50) 

In the case of an airfoil of nonzero thickness and arbitrary camber, the domain 
outside the airfoil in w = x + iy space will map onto a nearly circular domain in the 
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I 


z = £ + irj plane. The consequent introduction of elliptic coordinates will not help to 
eliminate the problem of mesh nonorthogonality, and one side of the computational 
rectangle will be an irregular line. 

The shape of the mesh cells at axial infinities (fig. 3) was determined explicitly 
(ref. 1) in order to make the application of the proper boundary conditions (eqs. (31), 
(39), and (44)) easier. 


Finite Area Method 

A uniform orthogonal mesh in the (X, Y) computational space remaps into a non- 
orthogonal mesh in the (x, y) physical space (fig. 3). The local, bilinear, isopara- 
metric mapping functions of the general form 


where 


b 4E v i+xx p>< i+y v 




Y = ±1 
P 


where p = 1, 2, 4, 5 


(51) 


(52) 


and b stands for x, y or G(x, y), will map each unit square mesh cell from the (X, Y) 
computational rectangle (fig. 5) into a corresponding distorted mesh cell in the (x, y) 
physical plane. 

The mass flux balance must be satisfied within each auxiliary control cell (ACC), 
which is represented as a shaded region in figure 5 and centered around each mesh 
point. This cell is formed from parts of the four neighboring elementary mesh cells 
(EMC). Each EMC was separately mapped from the (x, y) space to the (X, Y) space. 
Therefore the value of the desired parameter at the center of each mesh face of that 
ACC (points N, S, E, W) is evaluated as an arithmetic mean (ref. 3) of the four separate 
results that are calculated on the basis of the local mapping functions in each of the 
four neighboring EMC’s. 


Artificial Time Concept 

An attempt to solve the steady, full- potential equation - as an asymptotic solution 
of the exact, physically unsteady, full- potential equation for long times - would be 
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very uneconomical because of the small time steps required by the numerical stability 
considerations. Instead it is customary to use an artificial time- dependent equation 
(ref. 4) of the form 

( M r " 1 )^, ss ~ nn + 2o; 1 <P, st + 2a 2 < f > , nt + ecp , t = 0 < 53 > 

This time- dependent process does not model the true time evolution of the flow. 
The major advantage of this technique, however, is that the consistency and accuracy 
of the numerical method used in determining the transient solutions are irrelevant - 
the consistency being finally achieved asymptotically when the solution no longer de- 
pends on time. 

Artificially time- dependent differencing was used in a way suggested by Jameson 
(ref. 5), who considered iteration sweeps through the computational field as successive 
intervals in an artificial time. Using the successive line overrelaxation (SLOR) tech- 
nique requires introduction of a temporary or provisional value of the potential <p at 
every mesh point on the line along which SLOR is to be applied. The definition of such 
a term is 





(54) 


where u> is the overrelaxation factor (ref. 4), <pt . is the new value of the potential 

■*■>3 o 

(i. e. , one that will be obtained as the result of the current iteration sweep), and <p. . 
is the old value of the potential (obtained after the previous iteration sweep). Accord- 
ing to the numerical stability analysis (ref. 5), factor o> has values less then 2 when 
the flow is locally subsonic and equal to 2 when the flow is locally supersonic. 

The generation of the artificial time derivatives is illustrated in the following ex- 
ample. It is easy (refs. 4 and 1) to check that 





1 

(AX) 2 


represents 






(55) 


( 56 ) 
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and that the expression 



o o + , o \ 


1 

(AX) 2 


represents 



( 57 ) 


(58) 


where the superscripts E and H 
respectively. 

The vector of the corrections 


represent the central and upstream differencing, 


to the potential cp was determined from 
[«{Cj}={Rj} (59) 


where the tridiagonal [ 6] matrix of coefficients was determined from equations (22), 
(23), and (54). The vector of residues <R.> was determined from equations (21) and 

/rt ^ ' 


(24). 


COMPUTER PROGRAM 
General Description 

The computer program CAS2D consists of 10 routines and a separate input data 
operation (fig. 6). The input data are discussed in detail in the following section. 

The subroutine MAIN is the principal part of CAS2D in the sense that all other 
routines are called from that routine. MAIN reads the input data and rotates the air- 
foil to its actual stagger angle. This routine also calculates the length of the central 
slit in the circle-plane (fig. 4) and interpolates the symmetrically spaced points on the 
unwrapped surface of the airfoil in the (u v ) plane (fig. 4). Furthermore MAIN cal- 
culates (iteratively) the flow parameters at downstream infinity. It also determines 
the exact value of the circulation r (eq. (45)). Routine MAIN also tests the conver- 
gence rate of the iterative process by comparing the value of the expression (r n+ ^ 

- r n )/r n with the prespecified convergence- rate input parameter called CONVER. 
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Subroutine CONMAP iteratively performs the point- by-point conformal mapping 
(eq. (46)) from the (x, y) plane onto the circle-plane (£, 77 ) (fig. 4). Furthermore 
CONMAP "unwraps” the circle and calculates the elliptic polar coordinates u g and 
v e (fig. 4). 

Subroutine SPLIF fits a cubic spline through the lower boundary of the computa- 
tional domain in the (u v ) plane, which corresponds to the surface of the airfoil. 

g e 

Subroutine INTPL interpolates the values of the elliptic polar coordinates in the 
(u , v ) plane at points that are equidistantly spaced in the X-direction in the (X, Y) 
computational plane with respect to the image of the upstream infinity (X = 0; Y = 0). 
This is a necessary step in obtaining a grid that is periodic in the vertical direction in 
the (x, y) plane. 

Subroutine REMAP analytically determines coordinates of the mesh points in the 
physical plane; that is, REMAP performs a backtransformation process from the 
(X, Y) plane to the (x, y) plane. 

The mesh points defining the axial infinities should be positioned along the line 
x = Constant (fig. 3 and eq. (31)). Because of the way that the potential jump r is 
enforced across the cut, the points at x = +°° should be equidistantly positioned in the 
y direction with respect to that cut (fig. 3). 

Subroutine XYINF determines the coordinates of the axial infinity points explicitly 
because they cannot be obtained from the conformal mapping. XYENF also determines 
the coordinates of points in the imaginary rows and columns outside the actual compu- 
tational domain. 

Subroutine XSWEE P iteratively performs the actual flow calculation by using the 
SLOR technique. XSWEE P uses type- dependent rotated finite differencing (eq. (22)) 
and the artificial time concept in order to evaluate coefficients of the tridiagonal cor- 
rection to the potential matrix. The residues are evaluated by using the finite area 
technique (eqs. (20) or (21)). The artificial viscosity is added explicitly in conserva- 
tive form (eq. (24)). The form of finite differencing (i. e. , the fully conservative 
scheme (eq. ( 20 )) versus the quasi- conservative scheme (eq. ( 21 )) is indicated at the 
beginning of each of the two versions of the XSWEE P routine. The choice of the ver- 
sion of the XSWEE P routine that should be used depends on the user. Both versions 
give practically identical results for the moderately shocked flows. At the same time 
the fully conservative version takes about 50 percent more execution time. 

Subroutine BOUND applies boundary and periodicity conditions after each complete 
sweep through the flow field performed by XSWEE P. BOUND is also called after each 
mesh refinement. 

Subroutine CPMACH calculates a Mach number at every mesh point in the flow 
field and prints a Mach number chart after the iterative calculation process has con- 
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verged on the final grid. CPMACH also calculates and prints out values of the iter- 
atively obtained Mach number, density, and flow angle at the downstream infinity, 
which can be compared with the exact values obtained from MAIN. 

If the iterative calculation process is to be repeated on the next finer mesh, sub- 
routine MESH is called. MESH doubles the number of mesh cells in each direction 
(X and Y) and interpolates the values of the reduced potential G(x, y) (obtained on the 
previous grid) onto the new grid, thus creating an improved initial guess for the iter- 
ative process on that refined grid. 


Input 

The first card of the input data contains an arbitrary text with up to 80 characters 
describing the airfoil or the test case. The actual input parameters, which specify 
the flow and geometry, are given as real numbers. 

The second card contains XCELL, YCELL, and PMESH, where 

XCELL number of mesh cells on surface of airfoil for coarse grid. The number 
of mesh cells must be an even number; the suggested minimum value is 

XCELL = 16. 

YCELL number of elliptic layers of mesh cells enveloping airfoil, that is, half 
the number of mesh cells between two neighboring blades (for coarse 
grid). The suggested minimum value is 

YCELL = 4. 

PMESH total number of grids used. The present version of CAS2D is capable of 
calculating on a coarse grid and three consecutively refined grids. 
Hence the maximum allowable value is 

PMESH = 4. 

The third card contains ALPHA1, TWIST, and ALPHA2, where 

ALPHA1, angle, a_ 00 , (in degrees) between x-axis and free stream at up- and 

ALPHA2 downstream infinity, respectively (fig. 1) 

TWIST stagger angle, /3 (in degrees) between x-axis and airfoil chord line 

(fig. 1) 
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The fourth card contains PITCH, ROl, and R02, where 

PITCH gap- chord ratio, h/c (fig. 1). The iterative determination of certain 

parameters in the mesh- generating routines might fail for small values 
of PITCH because of the computer- dependent accuracy. Therefore it 
is suggested that CAS2D be used for cascades with 

PITCH > 0.65 

ROl, R02 radii (dimensionless) of airfoil at leading and trailing edges, respectively. 
For theoretically zero radii with respect to chord length, use 

roi = i. no -4 

R02 = 1. *10~ 4 

The fifth card contains FMACH, CONVER, and AR, where 

FMACH Mach number M of free stream at upstream infinity. The maximum 
allowed value is 


FMACH = 0.99 

This value should be lowered if it produces a shock that has a Mach 
number just ahead of it greater than ~1. 3. 

RLX over relaxation factor (eq. (54)) on coarse grid for case of locally sub- 

sonic flow. The suggested value is 

RLX = 1.75 

This value will be automatically increased by 4 percent on each con- 
secutively finer grid. In the locally supersonic flow, the relaxation 
factor is automatically taken as 

RLX = 2. 

AR ratio of inlet area to exit area of blade-to-blade passage (fig. 2), 

AR = (Hh^/hh^,) 
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The input parameter AR should provide the appropriate value of the 
downstream Mach number. For an airfoil with a closed trailing edge, 
use 


AR = 1. 

If the airfoil has an open trailing edge extending to downstream infinity 
(e. g. , shockless airfoils obtained from the design method of ref. 8), 
the approximate value of AR could be calculated from the known thick- 
ness S w of the trailing edge, 

AR = 1. + (S w / PITCH) 

The sixth input card contains TITR1, TITR2, and TITR3, where 

TITR1, maximum number of iterations on each of the first three grids, re- 

TITR2, spectively. These three parameters are important because they 

TITR3 serve as a convergence criterion in the case of nonlifting flows 

— fi 

(T < 1. *10 ). The suggested values that will provide results with 

engineering accuracy are 

TITR1 = 150. 

TITR2 = 60. 

TITR3 = 30. 

In the case that PMESH = 4, the maximum number of iterations on the 
fourth grid will be automatically taken as TITR3/2. 

The seventh input card contains POINTS, GAMMA, and CONVER, where 

POINTS number of input mesh points on surface of airfoil. For a nonsymmetric 
airfoil, POINTS must be an odd number (counting the trailing- edge 
point twice). For a symmetric airfoil only an even number of points 
defining its lower surface should be given as input (counting the trailing- 
edge and leading- edge points once only). The maximum total number 
of input points in the present version of CAS2D is 

POINTS = 129. 
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Note that the number of the input points on the airfoil lower surface 
does not have to be the same as the number of the input points on its 
upper surface. Before preparing the input data it might be helpful to 
consult the four examples shown in the appendix. 

GAMMA ratio of specific heats of working fluid 

CONVER circulation rate of convergence criteria in case of lifting flows 


CONVER = 


(r n+1 - T n ) 


-,n 


The suggested input value is 

CONVER = 1. *10" 5 


The x’ and y coordinates of input points on the surface of an airfoil in the cas- 
cade are given on the rest of the input cards starting with the eighth card. These co- 
ordinates will be normalized with respect to the airfoil chord length. The input co- 
ordinate system (x< , y ) could be arbitrarily positioned with respect to the airfoil 
(fig. 7). 

The input points are numbered in the clockwise direction starting from the 
trailing-edge point. The numbering must end with the same trailing-edge point for a 
nonsymmetric airfoil and with the leading- edge point for a symmetric airfoil. The 
coordinates are printed on input cards in such a way that x' coordinates of all the input 
points are given in the first column, followed by the corresponding y' coordinates in 
the second column (see appendix). 


Output 

The results of the CAS2D computer program appear in printed form only. There 
are no output files stored in the computer. 

First, the printout gives the name of the program, the name of the programmer, 
and the name of the airfoil or test case as specified on the first input data card. 

The four columns that appear next show normalized x» and y' coordinates 
prerotated to a zero stagger angle (in the first two columns) and as they appear 
after rotating the airfoil for the value of stagger angle TWIST (in the last two 
columns). 
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Next, the values of the flow parameters at the downstream infinity are given. 

They are obtained from an iterative procedure involving the input flow parameters 
and the mass conservation principle. 

It is possible to monitor the iterative process by using the following parameters, 
which are printed after each complete iterative sweep through the flow field. These 
parameters are 

ITER number of iteration sweep just completed 

IR, JR coordinates of point where residue had largest absolute value 

MAX* RESIDUE maximum residue (eq. (20) or (21) plus (24) in flow field. Its loca- 
tion is at the point (ER, JR). 

IC, JC coordinates of point where correction to potential had maximum ab- 

solute value 

MAX* CORRECT maximum value of calculated correction to potential. This correction 
was introduced at the point (IC, JC). 

CIRCULATION value of circulation r 

RELAX* COEF value of relaxation factor RLX used in last iteration sweep (eq. (54)). 

ISTG number of leading- edge stagnation point. From that point the next 

iteration sweep will start proceeding along the airfoil suction 
surface to the trailing edge and then again from ISTG along the 
pressure surface to the trailing edge. In such a way the problems 
of marching upstream in the locally supersonic flow and the con- 
sequent introduction of negative artificial viscosity are avoided. 

NSUP total number of supersonic points in flow field 

When the absolute value of the normalized convergence rate 

(r n+1 - r n )/r n 

becomes smaller than CONVER, the iterative process on that particular grid will 
terminate. For r = 0 (when the flow is nonlifting) the iterative process on each grid 
will terminate after ITER = ITRMAX on that particular grid. The printout will then be 
continued by listing the following values on the blade surface: 

X x 1 coordinate of point on airfoil surface. (The airfoil has been rotated by the 

stagger angle TWIST. ) 

Y y< coordinate of point on airfoil surface 
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XNORM x* coordinate normalized with actual chord 


CP coefficient of pressure 


C 


P 




DENS local density normalized with respect to density at upstream infinity 
MACH local value of Mach number 

Q/QINF ratio of local speed to speed at upstream infinity 


Q/QINF = [(u 2 



1/2 


If PMESH was given as different from PMESH = 1 in the input data, subroutine 
MESH will refine the basic mesh so that the new mesh will have twice as many mesh 
cells in the X and Y directions as the previous mesh. The printout will continue 
with a listing of ITER, IR, JR, . . . , etc. , on that new mesh. 

Finally, when the last iteration sweep is completed on the last specified mesh, 
the chart of the Mach numbers (multiplied by 10) in the entire flow field will be printed. 
The first horizontal line of numbers ranging from j = 2 to J = MAXY corresponds to 
the elliptic mesh layers enveloping the airfoil. Here J = 2 corresponds to the mesh 
layer along the periodic boundary; MAXY corresponds to the mesh layer along the 
surface of the airfoil. 

The first column of numbers on this chart designates the number of the mesh 
point, with 1=2 corresponding to the (lower) trailing- edge point. The Mach number 
chart will be deleted from the printout if the Mach number at upstream infinity is 

M_ oo < 0. 0011. 

The Mach number chart is followed by the values of the Mach number, density, and 
free- stream angle at downstream infinity calculated after the end of the last iteration 
sweep on the final grid. They are followed by the final results of the calculation, a 
list of the flow parameters on the airfoil surface and the C -distribution chart. 

r 
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RESULTS 


The computer program CAS2D was tested four times. All calculations were per- 
formed on three consecutively refined grids measuring 24 x 6, 48 x 12, and 96 x 24 
mesh cells, respectively. The relaxation parameter used on the coarse grid was 

RLX = 1. 65. The convergence criterion for the circulation rate was CONVER = 

— fi 

1. *10 on the coarse grid. The maximum allowable numbers of iterations on the 
three grids were specified as 

TITR1 = 120. TITR2 = 60. TITR3 = 30. 

The remaining input parameters for each test case are summarized in table L 
The last column in that table gives references that provide results widely accepted as 
being exact. A private communication from D. A. Caughey of Cornell University pro- 
vides the results for an isolated NACA 0012 airfoil in free air. Figure 8 shows that 
these results are in excellent agreement with the results obtained by using CAS2D for 
PITCH = 3. 6 (the nonlifting incompressible test, test case 1). 

Figure 9 gives the comparison between the second test case, the lifting incom- 
pressible test, and the results of reference 9. The agreement is again excellent des- 
pite the fact that the airfoil used is cusped. Figure 10 compares the nonlifting tran- 
sonic, shocked test case, test case 3, with the results obtained by D. A. Caughey. 
There is practically no difference between the results of the fully conservative and 
quasi- conservative schemes for shocks of that strength. This is in excellent agree- 
ment with the theoretical analysis performed by Caughey (ref. 10). 

The fourth test case was done on a shockless cascade of airfoils obtained by J. 
Sanz of Langley Research Center (ICASE) (private communication), who used a com- 
puter code developed by Bauer, Garabedian, and Korn (ref. 8). This airfoil had an 
open trailing edge extending to downstream infinity. After the trailing edge had been 
rounded (closed), the proper flow parameters at the downstream infinity were ob- 
tained by using AR = 1. 023 in the input. The agreement between the results of CAS2D 
and the results of Sanz are very good (fig. 11). The discrepancy of results near the 
trailing edge is a consequence of the geometric modifications in that region. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, March 5, 1980, 

505-32. 
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APPENDIX - LIMITATIONS OF CAS2D CODE 


Each user should be aware of certain limitations in the CAS2D code. A gap- chord 
ratio less then approximately PITCH = 0. 65 often leads to problems associated with 
the mesh generation routines. Keeping the stagger angle TWIST within the range 

+45° < TWIST < -45° 

I 

avoids large mesh distortions when PITCH is small. The Mach number at upstream 
infinity, FMACH, must be less than unity, and airfoil camber should not be excessive. 
Four example cases of input data are shown in figures 12 to 15. 

The format in which the input should be given is as follows: 


Card 

Format 

1 

2 to 7 
8 to end 

20A4 

9X, E12.5, 7X, E12. 5, 7X, E12. 5 
2F10.6 
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TABLE I. - INPUT DATA FOR FOUR TEST CASES 


Test 

case 

Airfoil 

Mach number 
of free stream 
at upstream 
infinity, 
FMACH 

Angle between 
free stream 
and x axis at 
upstream infinity, 
ALPHA1, 
deg 

Stagger 

angle, 

TWIST, 

deg 

Angle between 
free stream 
and x axis at 
downstream 
infinity, 
ALPHA2, 
deg 

Gap- chord 
ratio, 
PITCH, 

Ratio of 
blade 
passage 
exit area 
to inlet 
area, 
AR 

Comparison 

with 

reference 

1 

NACA 0012 

0.001 

0 

0 

0 

3.6 

1.0 

(a) 

2 

Gostelow 

(cusped) 

.001 

53.5 

-37.5 

30.0249 

.990157 

1.0 

9 

3 

NACA 0012 

.8 

0 

0 

0 

3.6 

1.0 

(a) 

4 

Jose Sanz 
(shockless) 

.711 

30.81 

-9.32968 

-.09 

1.02824 

1.023 

(b) 


a Private communication from D. A. Caughey of Sibley School of Mechanical and Aerospace Engineering, Cornell 
University. 

^Private communication from J. Sanz of ICASE, Langley Research Center. 



Figure 1. - Planar cascade of airfoils. 
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Figure 4. - Geometric transformation sequence, where TE denotes 
trailing edge and LE denotes leading edge. 
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Coefficient of pressure, 



Figure 8. - Comparison of non lifting incom- 
pressible test case (test case 1) with un- 
published data of D.A. Caughey of Sibley 
School of Mechanical and Aerospace Engi- 
neering, Cornell University. Airfoil, NACA 
0012; Mach number of free stream at upstream 
infinity, M^ (FMACH), 0.001; stagger angle 
between x axis and airfoil chord line, pfIWIST), 
0°; angle between free stream and x axis at 
upstream infinity, a_ 00 (ALPHAl), 0 P ; angle 
between free stream and x axis at downstream 
infinity, (ALPHA2), 0 P; ratio of inlet area 
to exit area of blade-to-blade passage, AR, 1. 0. 


O Test case 2 (finite area results}; 
-1.0— AR, 1. 0 

Results of ref. 9 



Figure 9. - Comparison of lifting uncompres- 
sible test case (test case 2) with results of 
referenced Airfoil, Gostelowcusped; Mach 
number of free stream at upstream infinity, 
M_oo (FMACH), 0.001; stagger angle between 
x axis and airfoil chord line, |3 (TWIST), -37.5°; 
angle between free stream and x axis at up- 
stream infinity, OLoolALPHAl), 53.5°; angle 
between free stream and x axis at downstream 
infinity, a +00 (ALPHA2), 30.0249°; gap-chord 
ratio, h/c (PITCH), 0.990157. 


to 

CD 






(a) Test case 3 (finite area results). (b) Unpublished data of D. A. Caughey. 

Figure 10. - Comparison of non lifting transonic-shocked test case (test case 3) with unpublished 
data of D. A. Caughey of Sibley School of Mechanical and Aerospace Engineering, Cornell 
University. Airfoil, NACA 0012; Mach number of free stream at upstream infinity, M_ oc (FMACH), 
0.8; stagger angle between x axis and airfoil chord line, |3 (TWIST), 0P; angle between free stream 
and x axis at upstream infinity, a -OQ (ALPHA1), OP; angle between x axis and airfoil chord line at 
downstream infinity, a +OQ (ALPHA2), (P ; gap-chord ratio, h/c (PITCH), 3.6; ratio of inlet area to 
exit area of blade-to-blade passage, AR, 1.0. 
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1.4 



x'/c 


Figure 11. - Comparison of lifting transonic 
shockless test case (test case 4) with unpub- 
lished data of J. Sanz of 1CASE, Langley 
Research Center. Airfoil, J. Sanz shockless; 
Mach number of free stream at upstream in- 
finity, M^IFMACH), 0.711; stagger angle 
between x axis and airfoil chord line, p (TWIST), 
-9.32968°; angle between free stream and x 
axis at upstream infinity, a^fALPHAD, 30.81°; 
angle between x axis and free stream at down- 
stream infinity, a +00 (ALPHA2), -0.09°; gap- 
chord ratio, h/c (PITCH), 1.02824. 


OOOO l 00 

MAC A 0012 

- INCOMPRESSIBLE 

NONLIFTING CASE 


0000200 

XCELL = 0 

. 20000D+02YCELL 

= 

0. 05000D+02PMESH = 

4. OOOOOD+OO 

0000300 

ALPHA 1= 0 

. OOOOOD+OOTW 1ST 

= 

0. OOOOOD+OOALPHA2= 

0. OOOOOD+OO 

0000400 

PITCH = 3 

. 6COOOD+OORO 1 

= 

0.01530D+OORQ2 

0. 00050D+00 

0000500 

FMACH = 0 . 00 1 OOD+OORL X 

= 

1.72000D+00AR 

1 . OOOOOD+OO 

0000600 

TITR1 = 1 

. 50000D+02TI TR2 

= 

0 . 60000D+02T I TR3 = 

0. 36000D+02 

0000700 

PO I NTS= 0 

. 18000D+02GAMMA 

■ 

1 . 40OOOD+OOC0NVER- 

1.00000D-06 

0000800 

1 . 000000 

0.000000 




0000900 

0.950000 

-0. 008070 




0001000 

0.900000 

-0.014480 




0001 100 

0. 300000 

-0.026230 




0001200 

0.700000 

-0.036640 




0001300 

0.600000 

-0.045630 




0001400 

0. 500000 

-0.052940 




0001500 

0. 400000 

-0.058030 




0001600 

0. 300000 

-0.060020 




000 1 700 

0. 250000 

-0.059410 




0001800 

0. 200000 

-0.057370 




0001900 

0. 150000 

-0.053450 




0002000 

0. 1 00000 

-0. 046330 




0002100 

0.075000 

-0. 042000 




0002200 

0.050000 

-0.035550 




0002300 

0.025000 

-0.026150 




0002400 

0.012500 

-0.018940 




0002500 

0.000000 

0.000000 





Figure 12. - Input data for incompressible non lifting test case. 



0000100 GGSTELOUI CUSPED AIRFOIL 
0000200 XCELL = O. 24000D+02YCELL 

0000300 ALPHA 1* 0. 53500D+02TWIST 

0000400 PITCH * 9.90157D-01R01 

0000500 FMACH = 0. OOlOOD+OORLX 

0000600 TITR1 = 1 . 50000D+02TITR2 

0000700 POINTS* 0. 39000D+02GAMMA 


0000800 

1.000000 

0.000000 

0000900 

0.993590 

0.000390 

0001000 

0.98S070 

0. 003090 

0001100 

0.918570 

0.017050 

0001200 

0.816750 

0.029510 

0001300 

0.691660 

0.035040 

0001400 

0.562250 

0.031740 

0001500 

0.425130 

0.020330 

0001600 

0.321980 

0.007730 

0001700 

0.272560 

0.001030 

0001800 

0.232030 

-0. 004630 

0001900 

0. 165360 

-0.013530 

0002000 

0. 110980 

-0.019350 

0002100 

0. 057580 

-0.021490 

0002200 

0.030890 

-0. 019220 

0002300 

0.015190 

-0.015100 

0002400 

0. 005990 

-0. 009930 

0002500 

0. 000010 

-0. 000040 

0002600 

0.001 130 

0. 008660 

0002700 

0.004610 

0.015330 

0002800 

0.017060 

0.023890 

0002900 

0. 035630 

0. 042360 

0003000 

0.059530 

0.055450 

0003100 

0.038440 

0. 067920 

0003200 

0. 122350 

0.079540 

0003300 

0. 161610 

0. 090090 

0003400 

0. 206980 

0.099280 

0003500 

0.259910 

0. 106710 

0003600 

0. 323170 

0. 1 1 177:0 

0003700 

0.403050 

0. 1 13030 

0003S00 

0.521140 

0. 106230 

0003900 

0.710520 

0.077590 

0004000 

0. 736060 

0.060730 

0004100 

0.826600 

0.050520 

0004200 

0.914960 

0.025760 

0004300 

0.961260 

0. 011740 

0004400 

0.9371 10 

0.003330 

0004500 

0.998680 

0.000330 

0004600 

1 . oooono 

0.000000 


= 0.06000D+02PMESH = 4.00000D+00 
=-0.37500D+02ALPHA2= 3.00249D+01 
= 0.O1230D+00R02 = 0.00005D+00 
= 1 . 75000D+00AR = 1. OOOOOD+OO 
= 0-50000D+02TITR3 = 0.30000D+02 
= 1 . 40OOOD+OOCONVER= 1 . 000000-05 


Figure 13. - Input data for incompressible lifting test case. 


0000 l oo 

NACA 001 

2 - TP ANCON I C NONLIFTING CASE 


0000200 

XCELL = 

2 4 0 0 0 D *■ 0 2 Y C E L L 

= 

0. 06000D+02PMESH = 

4. OOOOOD+OO 

0000 200 

ALPHA 1 = 

:> . OOOOOD+OOTW I ST 

= 

0. OOOOOD+OOALPHA2* 

0. OOOOOD+OO 

00» I04O0 

PITCH = 

3. 600000 +00 RQ 1 

= 

0.015S0D+00R02 = 

0. 00050 D+ 00 

0000500 

FMACH = 

3. 300000 +OORLX 

- 

1 . 73000D+00AR 

1 . OOOOOD+OO 

0000600 

TITR1 = 

1 . 50000D+02TI7R2 

= 

0. 60000D+02TITR3 = 

0. 36000D+02 

0000700 

POINTS* 

D. 13000D+02GAMMA 

= 

1 . 400C0D+00C ONVER= 

1 . 000 00 D- 06 

0000800 

1 . 000000 

0. 000000 




0000900 

0.950000 

-0.003070 




0001000 

0. 900000 

-0. 01 4430 




0001 100 

0 . 300000 

-0.026230 




0001200 

0. 700000 

-0.036640 




0001300 

0. 600000 

-0. 045630 




000 1 400 

0. 5000 00 

-0. 052940 




0001500 

0 . 400000 

-0.058030 




OOO 1 600 

O. 300000 

-0 . 060020 




000 1 *»00 

0. 2-50000 

-0.05 ->410 




000 t 300 

0 . 200* *00 

-0.057370 




000 1 900 

0. 150000 

-0. 05 3450 




0002000 

0 - 1 00000 

—0 . 04/. 8 30 




OiOO 2 1 00 

0. 075000 

-0 . 042000 




0002200 

0 . 050000 

-0. 035550 




0002300 

0. 025000 

-0.026150 




0002400 

0.012500 

-0. 013940 




0002500 

0.000000 

0. 000000 





Figure 14. - Input data for transonic, shocked non lifting test case. 



0000100 

JOSE SANZ 

SHOCKLESS CASCADE - TIP SECTION 




0000200 

XCELL = 0. 

, 24000D+02YCELL 

S3 

0.06000D+02PMESH = 3. OOOOOD+OO 

0006900 

-0. 007900 

0.006900 

0000300 

ALPHA 1= 0. 

, 30810D+02TWIST 

=- 

-9. 32968D+00ALPHA2=~0 . 09000D+00 

0007000 

-0.006400 

0.016500 

0000400 

PITCH = 1. 

. O44O0D + 0ORO1 

= 

0. 01000D+00R02 = 0.01000D+00 

0007100 

0.031 100 

0.055600 

0000500 

FMACH = 0. 

. 7 1 100D+00RLX 


1 . 68000D+OOAR = 1 . 026000+00 

0007200 

0.054800 

O. 072200 

0000600 

TITR1 = 1. 

■20000D+02TITR2 


0. 60000D+02T I TR3 * 0.30000D+02 

0007300 

0.082600 

0.089600 

0000700 

POINTS= 1. 

. 27000D+02GAMMA 

= 

1 . 4QOOOD+OOCONVER= 1 . OOOOOD-OS 

0007400 

0. 1 10300 

0. 105400 

0000300 

0.994000 

0. 171500 



0007500 

0. 136100 

0. 1 19000 

0000900 

0.990000 

0. 163000 



0007600 

0. 153600 

0. 130000 

0001000 

0.980000 

0. 160200 



0007700 

0. 173400 

0. 139100 

0001100 

0.955800 

0. 161200 



0007800 

0. 195600 

0. 146600 

0001200 

0.908400 

0. 163200 



0007900 

0.210900 

0. 152900 

0001300 

0 . 354600 

0. 165300 



0008000 

0.224500 

0. 153200 

0001400 

0.802700 

0. 166300 



0008100 

0.236800 

0. 162800 

0001500 

0.754500 

0. 167100 



0008200 

0.247800 

0. 166800 

0001600 

0.710900 

0. 166400 



0008-300 

0. 257900 

0. 170300 

0001700 

0. 671500 

0. 164800 



0008400 

0.. 267100 

0. 17.3300 

0001800 

0.635900 

0. 162600 



0008500 

0.275600 

0.176000 

000 L 900 

0.603500 

0. 159900 



0003600 

0.233400 

0. 173400 

0002000 

0. 573900 

0. 156800 



0003700 

0. 290700 

0. 130600 

0002100 

0.546700 

0. 153400 



0008800 

0.297500 

0. 132500 

0002200 

0.521700 

0. I 4*300 



0008900 

0. 303900 

0. 184300 

0002300 

0 . 493400 

0. 146100 



0009000 

0. 309900 

0. 135900 

0002400 

0 . 476800 

0. ) 42200 



0009100 

0.315700 

0. 137400 

0002500 

0.456600 

0 . 1 33 300 



0009200 

0. 321100 

0. 188300 

0002600 

0.437700 

0. 134400 



0009300 

0 . 326200 

0. 1 90000 

0002 700 

0. 419300 

0. 1 30400 



i>009400 

0. 331 100 

0. 191200 

0002300 

0.402900 

0. 126400 



000*500 

0. 335700 

0. 192200 

0002900 

0. 336900 

0. 1 22400 



0009600 

0. 340200 

0. [95200 

0003000 

0.372100 

0. 1 1310«‘» 



0009700 

0. 344500 

0. 194100 

0003100 

0. 357800 

0. 1 14200 



0009800 

0. 343600 

0. 195000 

0003200 

0. 343900 

0. 1 lOjuO 



0009900 

0. 352600 

0. 195300 

0003300 

0. 330700 

0. 1 O/ 400 



00 1 0000 

0. 356600 

0. 196500 

0008400 

0.313000 

0. 102600 



0010 100 

0 . 360500 

O. 1*7 300 

000 3500 

0. 305300 

0. 093800 



0010200 

0. 364400 

0. 1*3000 

0003600 

0. 294000 

0 . 095000 



0010300 

0. 368.300 

0. 1*3600 

0003700 

0. 282600 

0.091200 



00 1 0400 

0. 372200 

0. 19*300 

0003300 

0.271500 

0.087500 



0010500 

0. 376200 

o . 1 99900 

0003900 

0. 260700 

0 . 03 3800 



00 1 0600 

0. 380300 

0. 200500 

0004000 

0. 250200 

0.080200 



00 1 0700 

0. 398700 

0. 202700 

0004100 

0. 240100 

0. 076600 



'30 1 0300 

0. 408000 

0. 204000 

0004200 

0. 230100 

0.07 3000 



0010900 

O. 41 7500 

O. 204900 

0004 300 

0. 220400 

O.06"'500 



on i i ooo 

0 . 427700 

0 . 205800 

O0044O0 

0. 210700 

«■>. 066000 



00111 00 

0.433200 

0. 206700 

0004500 

0 . 20 1 400 

0. 062600 



00 1 1 200 

O. 449500 

0. 207400 

0004 600 

0. 191900 

0. 059300 



00 1 1 300 

0.461200 

0. 208000 

0004700 

0. 132700 

0. 055900 



O01 1400 

0. 473700 

0. 200 500 

0004800 

0. 17 3600 

0.052500 



00 1 1 500 

0. 486800 

0 . 208800 

0004900 

0. 164700 

0.049100 



00 1 1 600 

0. 500800 

0 . 209 1 00 

0005000 

0. 155700 

0.045700 



001 1700 

0.515600 

0. 20 71 200 

0005 1 00 

0. 1 46900 

0. 04240' > 



0011800 

0. 5.31500 

0. 209100 

0005200 

0. 1 38000 

0. 03*000 



001 1900 

0. 548600 

0 . 208800 

0005300 

0. 129200 

0. 035800 



00 1 2000 

0.5671 00 

0. 208200 

0005400 

0. 120 300 

0. 032400 



OO 1 2 1 00 

0. 537200 

0. 207400 

(-1005500 

0.11 1600 

0. 029 200 



0012200 

0 . 609400 

0. 206 400 

0005600 

0. 102700 

0.025900 



0012300 

0. 6.33900 

0 . 205000 

0005700 

0.094000 

0.0227O0 



0012400 

0.661500 

0 . 20 3200 

0005300 

0.085100 

0.019400 



00 1 2500 

0.692500 

0. 201000 

0005900 

0. 076400 

0. 016400 



00 1 2600 

0. 727900 

0. 193500 

0006000 

0. 067500 

0.013 300 



0012700 

0.763400 

0. 195500 

0006100 

0. 058900 

0. 010400 



00 1 2800 

0.815100 

0. 1*2300 

0006200 

0.050000 

0 . 007700 
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Figure 15. - Input data for transonic, shockless lifting test case. 
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